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Abstract 


The goal of this investigation is to numerically simulate the effects of oscillatory 
actuators placed on the leading edge of an airfoil, and to quantify the effects of oscillatory 
blowing on an airfoil stall behavior. It has been demonstrated experimentally that periodic 
blowing can delay flow separation at high angle of attack. The computations are to be performed 
for a TAU 0015 airfoil at a high Reynolds number of -IxlO 6 with turbulent flow conditions. 
The two-equation Wilcox k — co turbulence model has been shown to provide reliable 
descriptions of transition and turbulence at high Reynolds numbers. The results are to be 
compared to Seifert’s experimental data. 


Introduction 

Studies indicate that when periodic flow excitation is applied to airfoils, it can effectively 
delay flow separation. There have been several approaches introduced to reduce flow separation, 
namely, deformation of the leading edge of an airfoil, steady blowing and several different 
methods of oscillatory blowing applied to the airfoil. This research is focused on steady or 
oscillatory blowing on the leading edge of an airfoil. Studies have indicated that 
oscillatory/periodic blowing is significantly more efficient than steady suction or steady blowing 
at controlling separation. 11 ,2,3] Not only that, the implementation of a controlled periodic 
excitation is more feasible than steady blowing. 

Experiments on a NACA 0015 airfoil with low Reynolds number have shown that the 
incompressible maximum lift coefficient, Ci, max , can be increased by more than 15%. In addition, 
the post-stall lift can be increased by as much as 50% and post-stall drag can be reduced by more 
than 50%. [41 The airfoil that is currently being investigated is the TAU 0015, which was tested in 
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a low-speed wind tunnel at Tel-Aviv University (TAU). Finally, experiments performed over a 
wide range of Reynolds and Mach numbers showed that periodic excitation can effectively delay 
turbulent boundary layer separation and reattach flows. [U,3] Further research in active control 
has shown that Reynolds numbers corresponding to a jetliner at takeoff condition were 
successfully applied.^ The numerical model being used is the compressible Navier-Stokes 
OSC2D code implementing the two-equation Wilcox k-a) turbulence model. Predictions will be 
compared with the experimental results produced by Seifert [5] as well as with the computational 
results from Joslin et al J 61 and Donovan et al. . 


Airfoil Geometry 

A common practice in experiments is to use a blunt or finite thickness trailing edge. 
However, for computational and grid generation purposes we implement a sharp trailing edge. 
The NACA OOxx airfoils with a finite thickness trailing edge are described by 
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where j/ is the thickness-chord ratio (e.g., l / = 0.15 for a NACA 0015 airfoil). The constant 

coefficients for the equation (1) are: ao — 0.2969, aj = -0.1260, aj — -0.3516, a? = 0.2843 and 
a 4 = -0.1015. [S] The coefficient of the x 2 component is changed to a 2 = -0.3537, to convert the 
NACA OOxx airfoil from a blunt trailing edge to a sharp trailing edge. 

The coordinate data for the TAU 0015 airfoil provided by Seifert includes a notch at the 

leading edge, a notch at x / = 0.76 on the upper surface, and a finite thickness trailing edge. The 
notch on the upper surface will be ignored in this study. By using a similar approach to that 
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described above, the TAU 0015 data provided by Seifert was converted to a sharp trailing edge 
by the transformations 


y-y^-o.ooisvsix, f 
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x + a00328 '\ 
1.00328 


is a function of the specified length for the upper surface of the 
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where x, = 


x + 0.00656 


is a function of the specified length for the lower surface of the 


1.00656 

\ / 

airfoil. The difference between equations (2) and (3) is due to the definition of points on the 
trailing edge at the upper and lower surfaces. 

Now, by modifying the NACA 0015 equation we can closely match the values produced 

from the equations above with y = y and x = x / , and let Co = 1.0 for comparison purposes. 

/ c o / c o 


First, the upper surface of the airfoil is described by 

y lt = 0. 75{a 0 x^ 2 + a t x + a 2 x 2 + a 3 x 3 + a 4 x 4 ), for 0.3<x<1.0 


and 
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. The lower surface is described by 
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where x, = (r + 0.00656] : . The coefficients a, are determined by least square fits of 

2 v 1 0.30656 


equations (4) to the respective portions of the sharp trailing edge TAU 0015. 

The relative differences in y between the TAU 0015 sharp trailing edge airfoil described 
by equations (1 )-(3) and the curve-fitted TAU 0015 airfoil described by equations (4) are less 
than ±0.5% with the exception of a couple points located near the discontinuity located at the 
notch on the leading edge, where the absolute difference in y were less than 0.0001. 

Now, we normalize the current curve-fitted TAU 0015 airfoil for computational purposes. 
The x-transformations are needed to make the leading edge (of the lower surface) coincide with 
the origin. We do this normalization and transformation by 


. _ x + 0.00656 
X L00656 


(5. a) 


for both the lower and upper surfaces, then 

y, 

* = 1 1.00656 (5.b) 

y_k 

. 1.00656 

again, for the lower and upper surfaces which was defined in equations (4). The purpose of using 
a curve-fitted equation is to define the surface boundary of the airfoil continuously with the 
exception of the notch at the leading edge. 
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Equations of Motion 

For compressible turbulent fluid flows, the equations of motion and k-Q) turbulence 
model equations written in terms of Favre mass-averaged quantities are 


Mass conservation: 


«£+ 9_(p, ),o 

dt dXj ^ j ’ 


Momentum balance: 
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Mean energy conservation: 




Turbulent mixing energy: 
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Specific dissipation rate: 
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where t is time, x, is the position vector, m, is the velocity vector, p is density, p is pressure, p is 
molecular viscosity, f y is the sum of molecular and Reynolds stress tensors, and qj is the sum of 

the molecular and turbulent heat flux vectors. In Equation (8), the quantities E = e + k + u,u> / 2 
and H =h + k + u,u i /2 are total energy and enthalpy, respectively, where h = e + p / p, and e and 


h denote internal energy and enthalpy, respectively. Additionally, is the Reynolds stress tensor. 
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The turbulent mixing energy k and the specific dissipation rate a) are needed to define the 


eddy viscosity 


H T =a 


[9] 


pk_ 

CO 


( 11 ) 


The total viscous stress tensor is 
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where the mean strain-rate tensor S, y is 
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We invoke the Boussinesq approximation that the Reynolds stress tensor is proportional to the 
mean strain-rate tensor 
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Finally, the heat flux vector q, is approximated as 

^ dh 


<li=- 


P + Pt 


dXj 


Pr L Pr T 

where Pr L and Pr, are the laminar and turbulent Prandtl numbers, respectively. The closure 
coefficients 191 are (5 = 3/ 40, (3* = 9 / 100, a = 5 / 9, a = 1, <y = 1 / 2 and c =1/2. 


(15) 
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Governing Equations 


The method used to formulate the flow governing equations past an airfoil makes use of 
an inertial reference frame in which any movement of the airfoil is simulated by multiple grid 
generation as the airfoil moves. The way to get around the multiple grid generation if the wing 
movement is predefined and the body is rigid. In this case, a single grid is generated for the 
problem and new grids are obtained by solving a set of equations based on rigid body rotation as 
the airfoil goes through its pitching motion. However, if the body is flexible or deforming, this 
method will not work. Hence, the multiple grid generation used to track the wing motion is more 
general and can handle both flexible and deforming bodies. The latter approach, combined with 
a moving grid, has been adopted here. 

Arbitrary shaped bodies present a difficulty for finite difference formulation of the 
problem in Cartesian coordinates. So, we need to map the physical space into a uniform 
computational space using a general curvilinear coordinate transformation. For the problem 
under study, the grid is moving in the physical space, so time must be included as one of the 
variables in the general coordinate transformation. 

Since wish to simulate deep dynamic stall, the full Navier-Stokes equations in 
nondimensional and conservative-law form are written as 

a + £, + ?;= Re-' (£„, + *•„) (16) 
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where Re is the Reynolds number, Pr is a Prandtl number and y is the ratio of specific heats (1.4 
for air). The speed of sound a is given for ideal fluids by a 2 = yp/p . The dynamic viscosity p 
is calculated from Sutherland’s formula and Stokes hypothesis A + (2/3) p = 0 is assumed. 


Pressure is related to the conservative flow variables, Q , by the equation of state 


P = {Y~l) 




(17) 


Equation (16) is based on an Eulerian description of motion. In this approach, a fixed 
control volume is used and the changes to the fluid are recorded as the fluid passes through the 
control volume. The choice of nondimensional parameters is arbitrary. The variables 
p (density), u andv (Cartesian velocities), e (total energy), t (time), and p , are scaled by 
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where °° refers to free stream quantities and l represents a reference length such as the chord of 
an airfoil. Note that the Reynolds number uses freestream speed of sound a„ rather than the 
usual , therefore it must be scaled by . For the remainder of this development the 

~ will be dropped for simplicity. 


Curvilinear Coordinate Transformation 

The Navier-Stokes equations can be transformed from Cartesian coordinates to general 
curvilinear coordinates where 
T = t, 

§=§(*,*/), ( 18 ) 

ri=r](x,y,t). 

The coordinate transformation, or computational space, introduced here are chosen so that the 
grid spacing in the curvilinear space is uniform and of unit length. This produces a 
computational space £ and t\ which is a rectangular domain of uniform unit cells. Therefore, 
the standard unweighted differencing schemes can be used in the numerical formulation. There 
will be a one to one correspondence between a point in computational space and a point in 
physical space except for regions where there are cuts due to topology or singularities. In those 
cases, it may be necessary to map one physical point to many computational points that normally 
occurs at computational boundaries. 
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The mapping relations are constructed through chain rule expansions, which are used to 
represent the Cartesian derivatives d x and d y of Eq. (16) in terms of the curvilinear derivatives 

where in matrix form 

A1 [i i n,l|A 

d, = o L n x 3^ ( 19 > 

aJ l° ^ 

In most cases the transformation from physical space to computational space is not 
known analytically, rather it is generated numerically. This is accomplished by reversing the role 
of the independent variables in the chain rule formulas, Eq. (19), to get 

a r =a,+* r a,+y T d y , 

3{ ~ X i^x + y^y’ ( 20 ) 

=^A+)’A 

and in matrix form 

a, I 1 3^r 9, 

0 *1 y ( d* W 

aJl° yJlA- 

Solving Eq. (21) to expand curvilinear derivatives d T , d 4 and d n in the terms of Cartesian 
derivatives can be written in matrix form 

3,1 [l 3V d T 

d x = 0 x f y 5 d { (22) 

_a v J [o x n y n \ [a„_ 

By evaluating equations (19) and (22) we get 

S, =-x£,-y£ y 
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V,=-x t r} x -y T ri y 

L = Jy n (23) 

r i x = -Jy i 

Vy = 

where the metric Jacobian, J , is defined by 

J~' = x i y„ -*„y 4 (24) 

and x r , y T are grid velocity components. Equation (16) after coordinate transformation using 
chain rule becomes 

Qr +iQi +n,Q n +iA + 1 h F n + %A +r i> F > i 

r n (25) 

= Re-'[«,(£,) { +n,(E.), +{,(*;), +<I,(F.),] 


Invariants of Transformation 

At this point of the development of the governing equations one can stop and use Eq. (25) 
to solve the flow field. It should be noted that Eq. (25) is not in a full divergence form, since 
metrics appear as coefficients in the above differential equations. Therefore, the flow properties 
are not strongly conserved. There are researchers who advocate the use of so-called “chain rule 
form” in problems such as shock capturing. Here, we shall restrict ourselves to the strong 
conservation law form, which will be derived below. Let us produce the strong conservation law 
form using Eq. (25). First we multiply this equation by 7~' and use the chain rule on all the 
terms, such as 
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If the latter part, Eq. (28), is eliminated then the strong conservation law form of the 
equations results in Part 1 equaling zero. For example, take the following term in Eq. (28) 


r ZA ,(n' 


\ J )i 
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and substitute for £ and r\ x to get 


(29) 

Now analytically differentiation is commutative and each of the above terms in then sums to 
zero, thus eliminating Part 2 or Eq. (28). Then equation Eq. (27) represents the strong 
conservation law form of Eq. (25). Further simplification of Eq. (27) using U and V , 
contravariant velocities, results in 
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with the contravariant velocities defined as 


U=£,+Z x u + Z y v, 

V =T! t +Tl x U + 11 y V 


and the stress terms are 


= | P [2 (&« { + r] x u n ) - (Zm + T] y v n )] , 
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Equation 30 represents the transformed governing flow equations in their final form. The 
implementation of the numerical method and its solution procedure are presented below. 
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Boundary Conditions 


Probably the most important aspect in successful application of any numerical technique 
is the proper treatment of boundary conditions, including impermeable boundaries such as solid 
wall and permeable boundaries such as inflow and outflow boundaries, which allow mass across 
their surfaces. 

A slip boundary condition is assigned on each tunnel wall, since the boundary layers on 
these walls cannot be resolved due to the limited number of grid points. In computational space, 
it is relatively easy to compute the velocity on the surface by means of extrapolating the 
contravariant velocity components. This procedure not only extrapolates the velocity 
components but also extrapolates the associated local contravariant base vectors which results in 
severe mesh-dependency of the solution. To correct the problem, it is not necessary to 
extrapolate the contravariant velocity but rather the tangential velocity components (physical 
space). This is done by linear first-order extrapolation, which sets the normal gradient of the 
velocity and other flow parameters equal to zero. For example, on an rj -constant plane the 
velocity component u is computed by setting 

— = Vw n = 0 (31) 

dn 

where n is a unit normal vector defined as 


yflf+vl 


(32) 


Combining Eq. (32) into Eq. (31) we get 


(33) 
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Equation (33) is solved for u implicitly in the £ -direction. Similarly, other variables can 

be extrapolated. On the surface of the wing, the no-slip boundary condition is specified. From 

Eq. (30) it is apparent that we have four dependent variables at every grid point. In order to 

update these variables at the surface, four equations are required. To obtain the correct velocity 

of the surface (wing), the contravariant velocities 

U=Z x {u-x T )+Z y {v-y 0 (34) 

V=rj x (u-x T )+ri y (v-y r ) 

are set to zero at the body (rj -constant plane). To extract the Cartesian velocity components w 
and v from Eq. (34), it is necessary to do a backward transformation from computational space 
to physical space, using 



% 

~v x L 



(35) 


A relation for pressure along the body surface is obtained by combining the two 
transformed momentum equations 

K (v 2 x + rfy f = (S A + ZyVy ) Pi + + P'y ) Pn (36) 

= p[d T n,+ud T r] x +vd t r} y ] 

The above relation for pressure boundary condition was found to produce virtually the same 
results as the simpler dp/dn = 0 condition in this study. However, relation (36) is more stable. 
Equations (34) and (36) provide a total of three equations. The fourth equation is obtained 
through specifying either constant wall temperature or an adiabatic wall. For the latter case 
which is used here, 
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must be satisfied at the wall. Assuming a perfect gas law, we have the following relation 
between pressure, density and temperature 
p = pRT . 

Differentiating the above relation with respect to ti and using the adiabatic condition and the 
assumption that dp/dn is zero between the wall and the first line of discretization adjacent to it, 
we have 



This equation can be written in computational space as 

(£ A + ) Pi + (ri +r ll)Pn =0 < 37) 

To obtain pressure and density. Equations (36) and (37) are implicitly solved in the ^ -direction. 
In the case of constant surface temperature, density at the wall is computed using the equation of 
state as 


where subscript s denotes properties at the surface. Now, with the addition of Eq. (37), there are 
five equations. Therefore, all conservative variables can be computed at the wall. 

For the inflow boundary, conservative variables are set to free-stream conditions. First-order 
extrapolation is used to update variables on the outflow boundary. 

There is a fictitious boundary condition or branch cut introduced by the coordinate 
transformation. The branch cut, which in this case is a plane, as shown in Figure 1, starts from 
trailing edge of the wing and extends up to the downstream boundary. Every point on this plane 
is represented by two points in the computational space. To correctly treat this boundary, it is 
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necessary to apply a periodic condition at the branch cut. This would result in a periodic block 
tridiagonal matrix which needs to be inverted at each time step. Periodic block tridiagonal 
inversion is much more costly than standard block tridiagonal inversion. Also, accuracy would 
require grid smoothness across the branch cut, which is difficult to obtain. However, due to the 
nature of the grid there is sufficient point clustering near the plane making it possible to 
extrapolate properties from interior nodes to two corresponding points in the computational 
space. Since every two points represents one point in a physical space, values there must be the 
same. Thus, an averaged value between the two corresponding points is computed and assigned 
to both points. 

The boundary conditions described in this section have not been incorporated into the 
implicit formulation, but are lagged behind one time step. The resulting code is flexible and can 
be adapted to different problems with a moderate amount of work. However, it is difficult to 
judge what happens to the practical accuracy of the method at the boundary due to this type of 
implementation. Visbal & Shang [l01 have evaluated implicit Navier-Stokes solvers in two- 
dimensions. They showed that for unsteady flow past a cylinder both implicit and explicit 
implementations of boundary conditions over the range of At (0.005-0.01) produced virtually 
identical results. 

For inviscid calculations, the slip condition is imposed on the surface of the wing. This is 
done in the computational space with the use of contravariant velocity components. Pressure is 
implicitly evaluated on the body using the unsteady form of the inviscid momentum equation. 

(v;+ri 2 y fp n = 

p[(v,l+ u { n x X + V (rj y ) t ]-pu (v ? - + Vs ) - P v (w + % u n ) = 
(Zjl x +Z,ri,)Pt+(tf+VI 2 ,)Pi 
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Furthermore, the pressure on the no-slip portion of the slot at the leading edge is 

(£+{}fp,= 

p[(S,), +t,'’ t )-pv(U + «a)= 

(i: +i:)i‘ t +{L’h+z,ri.)n, 

Free stream stagnation enthaply is held constant along the body. Using the equation of 
enthalpy. For the far-field, the non-reflecting characteristic boundary condition is used. This 
boundary condition requires a mapping between the normal and tangential component of 
velocities and the Cartesian component of velocities. Dr. Steinberg used the Symbolic 
Manipulation code Macsyma to construct this mapping 1 " 1 . 

For viscous calculations, the no-slip condition is applied on the surface of the wing by 
setting the contravariant velocity components to zero. Pressure is evaluated in a same manner as 
in the inviscid case, using the inviscid momentum equation. For the wall temperature, either an 
adiabatic condition or constant wall temperature are used. Density is computed using a zero 
gradient condition normal to the surface for an evaluated using zeroth-order extrapolation and for 
the far-field a non-reflecting characteristic boundary condition is used. The wake is evaluated by 
using an overlapping scheme for both the inviscid and viscous cases. 

Grid Topology 

A two dimensional C-grid will be used in this investigation of the TAU 0015 airfoil due 
to the high Reynolds number. The C-grid by its nature is not suitable for periodic boundary 
conditions in the ^-direction and to apply a true periodicity condition across the branch cut in the 
C-grid, the t^-sweep of the approximate factorization algorithm must be modified. The 
numerical finite difference code OSC2D implements the block pentadiagonal inversion and the 
cost is approximately less than two times the standard approach. The solution to the periodicity 
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problem developed by Visbal and Shang [12] is to overlap the branch cut region in the 
computational domain. This approach was verified in solving a 2-D Poisson equation in fully 
curvilinear coordinates using both O and C type grids. The results have shown that the concept 
is sound and can be applied to airfoil computations using the C-grid. Based, on the 
experimentation, it was determined that a five point overlap in the wake is sufficient to produce a 
smooth solution across the cut. [l 1] For the TAU 0015 airfoil, the grid topology is shown in 
Figure 1. The overlapping process in Figure 1 is shown by letting point B be the same point as 
B’, point A is A’ and so on. The polygon region of B-C-D-E-F-G-B is the iblank region of the 
grid. The term iblank indicates that a specified area of the grid is voided and no calculations are 
performed in that region. 


Approach and Discussion 

The experiments conducted by Seifert [5] used actuators at several locations operated at 
various frequencies with Mach number 0.15 and Reynolds number 1.2 million. However, this 
report is focused on the TAU 0015 with the actuator located at the leading edge with different 
frequencies and blowing coefficients. 

The different actuator boundary conditions that are being considered are introduced by 
the following oscillatory forcing functions 


m= 



sin 2 (n^) 


(15) 


where £ is the tangential direction to the surface. The jet suction/blowing is introduced through 
the velocity at the surface. 
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(16) 


M^ =o ’0= Vs +V%) sin ( F+ ^^) fit)** 

where r\ denotes the direction normal to the surface and e Jtl determines the angle the jet makes 
with the surface. The spatial variation of the jet is specified by /(£) and determines an effective 
jet width, 

H = \[f{^)] 2 d^, (17) 

where the integral is evaluated across the jet (0 < £ < a). The non-dimensional frequency is 
F + - fc/U„ , which relates the period of the jet cycle to the convection time of the flow over the 
airfoil. The steady momentum blowing coefficient, c M = 2(H/c^J jel jU „ J , and the oscillatory 

momentum blowing coefficient, (c M \ = 2{H / c^u' Jet J . determine the amplitude of the jet, 

where U jel is the mean jet velocity and is the RMS velocity of the jet oscillations. The 

flow profile that behaves better and produces a nearly identical flow [7] is sin 2 (n%). The rest of 
the boundary on the airfoil surface is prescribed as a no-slip boundary. 

In a previous study, predictions with an unstructured coarse grid showed some 
disagreement with measurements in the highly separated region in stall. [61 Furthermore, 
computations performed by Donovan et al J 61 and Joslin et al. [7] had trouble in predicting AC, 

versus where the experimental relation is monotonically increasing while the the 

numerical predictions are not. The approach herein is to use a fine, structured C-grid with the 
OSC2D finite difference code. The discontinuity of the blowing notch at the leading edge has 
caused some problems that will be overcome using the commercial package Gridgen to generate 
the geometric topology of the airfoil. The grid generated by Gridgen near the notch in the 
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leading edge is shown in Figure 3. The OSC2D code applies the two-equation Wilcox k-co 
turbulence model, with the assumption of a turbulent boundary layer starting from the leading 
edge. Furthermore, the two-equation k-co model has been shown to provide reasonably 
accurate description of turbulent and transitional flows in adverse pressure gradients. 19,131 
Computational results including the lift and pressure coefficients versus the normalized chord 
position and the angle of attack will be compared with Seifert’s experiment, as well as with the 
predictions of Donovan et al. [6] and Joslin et al. t?1 

For all computations in this study, the equations of motion have been integrated through 
the viscous sublayer. The Wilcox k-co model requires no special viscous damping functions to 
make possible integration through the sublayer. Viscous corrections by Wilcox 191 are needed 
only to improve transition predictions and to simulate subtle features of sublayer structure. The 
mean-flow properties are set according to no-slip and adiabatic-wall conditions. The no-slip 
condition also dictates that k vanishes at solid boundaries. Furthermore, we have implemented 
the slightly-rough surface boundary condition for to . Finally, the blowing portion of the slot is 
treated differently, where k and co are evaluated using zeroth-order extrapolation. 

Results are presented for oscillatory blowing at 4° angle of attack. Other angles of attack 
are being computed and will be published along with a comparison of the increase in lift 
coefficient provided by the oscillatory blowing. 

Figures 2a-e show the oscillatory blowing on the leading edge at the different times in 
one cycle beginning with zero blowing velocity, progressing to maximum blowing at one-quarter 
cycle, back to zero blowing velocity at one-half cycle, then to maximum blowing at three-quarter 
cycle, and finally back to zero blowing velocity to complete the cycle. This cycle begins after 
one full cycle has been completed and the results have reached periodic conditions. Fig. 2a 
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shows separation off the leading edge of the slot with a large recirculating zone. Fig. 2b shows 
strong vortex development due to shear between the flow over the leading edge and the blowing 
from the slot. A larger view shows these vortices being convected downstream where they have 
significant influence on separation on the upper surface of the airfoil at higher angles of attack. 
Fig. 2d shows strong interaction between the flow over the leading edge and the region above the 
slot where fluid is being removed by suction. 
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Figure 1. The coordinate transform between physical domain and computational domain for 
TAU 0015 airfoil. 
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Figure 2a. Velocity vectors at the leading edge of the TAU 0015 airfoil at the beginning 
of an oscillatory blowing cycle (zero velocity at slot). 



Figure 2b. Velocity vectors after one-quarter cycle of oscillatory blowing (maximum 
blowing at slot). 
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Figure 2c. Velocity vectors after half cycle of oscillatory blowing (zero velocity at slot). 



Figure 2d. Velocity vectors after three-quarter cycle of oscillatory blowing (maximum 
suction). 
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Figure 2e. Velocity vectors after a full cycle of blowing (zero velocity at slot) 



Figure 3. The grid structure at the leading edge of the airfoil. 
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